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ABSTRACT 

The earth tidal tilt data consists of two parts: one is 

the theoretical part caused by the tidal potential of the 
moon and the sun; another one is the seismic noise part that 
originates from meteorological, hydrodynamic, and cultural 
sources. The theoretical part is completely determinable 
by astronomical data. The seismic noise may be considered 
as a stationary stochastic process over short periods of 
time when there have been no great changes in the source 
mechanisms. Based on the above assumption, prediction and 
interpolation of the earth tidal tilt data can be achieved 
by solving the Wiener-Hopf equation. A computer program 
"PREDICT" is derived by using Levinson's recursive method to 
solve the discrete-time Wiener-Hopf equation. Test results 
on real data are satisfactory. Besides the earth tidal tilt 
data, the "PREDICT" program can also be applied to any 
stationary stochastic discrete-time signal. 
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INTRODUCTION 



I . 

The essential problem of earthquake prediction is to 
identify and monitor those physical characteristics of the 
failure process that uniquely herald the magnitude, location 
and origin time of an earthquake. Identification of these 
characteristics has tended to be heuristic and statistical 
in approach; however, evidence is steadily accumulating 
toward the support of a sound physical basis to earthquake 
prediction. Models of this failure process are becoming 
sufficiently detailed to effectively define the spatial and 
temporal character in each of the variety of geophysical 
measurements conducted in the epicentral region. 

This paper contains the theory and applications of a 
single-channel predictor for tilt at tidal sensitivity. 

Tidal tilt is caused by the direct tidal forces of moon and 
sun on the earth. There are also other sources such as 
meteorological activity and the ocean-loading effect which 
can cause significant tilts exceeding the magnitude of the 
primary source. The primary sources are predictable and, 
to a lesser degree, the effects of ocean and meteorological 
activity are also predictable. Other factors such as effects 
of local earthquake activity are random in nature and cannot 
be determined in advance. Therefore, we divide tidal-tilt 
information into two classes: one is predictable, the other 

random. In the nomenclature of communication theory, we 
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call the former signals and the latter noises. However in 
subsequent discussion we will, in fact, use the definition 
that the signal we seek is random in nature and the noise 
is the periodic but undesired part of the time series, 

A single-channel predictor in time series analysis uses 
past data to predict future data. The formulation of a self- 
predictor is basically one of designing a filter whose output 
in the time domain is a positive time-shifted version of the 
input, The prediction error is a measure of the stationarity 
of the time series. If the time series spectral density of 
equal lengths of data segments is a constant, then a filter 
designed from any segment of past data can also be used to 
predict future data from any other data segment. For the 
geophysical problem at hand, we derive such a filter during 
earthquake-free periods, because we assume that the failure 
process is, mathematically speaking, associated with non- 
stationarity in the data. The strategy is thus one of 
deriving a filter free of this "signal"; using that filter 
for subsequent incoming data, and monitoring the magnitude 
of the prediction error as a means of identifying onset of 
the failure process. 
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II. THEORY OF SINGLE-CHANNEL PREDICTOR 



A. BASIC TERMINOLOGIES 
1 . Digital Signals 

In order to use high-speed digital computers to 
analyze seismic signals, the first thing we must do is to 
convert the analog signals to digital signals. This involves 
the work to extract the values of the analog signal at 
discrete time points which are usually equally spaced. It 
is convenient to choose the time unit equal to the time- 
spacing At between two consective discrete time points so 
that the discrete time variable t takes on only integer 
values. The numerical values of the signal together with 
their associated discrete time points make up a digital 
signal. One example of a digital signal is 

Signal: 1, 3, 2, 4, 3, 5, 2, 4, 

Time t: •*♦, -2, -1, 0, 1, 2, 3, 4, 5, 

Because we have required time t to assume consecutive integer 
values, we can write the above digital signal simply as 

x = (•••, 1, 3, 2, 4, 3, 5, 2, 4, •••), 

t 

where the arrow is used to denote the time origin t = 0 and 
x represents the entire signal. Under this notation, x^ 



10 



will designate the amplitude of the digital signal x at 
time t. Thus in the above signal x^ = 2, x^ = 4 , and x ^ = 3 . 

2 , Wavelets 

Wavelets are a special class of digital signals. A 
wavelet is characterized by two properties, namely, 

(a) The one-sided property: A wavelet has a definite 

origin (or arrival) time in the sense that all values 
of the wavelet before its origin time are zero, that 
is b t = 0 for t < 0. 

(b) The stability property: A wavelet has finite energy 

so that it is a transient or dying-out phenomenon, 
that is 

oo 

£ b ? 

t=o z 

is finite. 

From the above two properties we know that a wavelet 
is a one-sided energy signal, It may damp out completely, 
i.e., become zero, after a certain time. This is the so- 
called finite-length wavelet. One example is 
C", 0, 0, 0, 5, 2, 3, 0, 0, 0, ••*) 

where we can simply write as 
(5, 2, 3). 

We call it a three-length wavelet. A wavelet can also be of 
infinite-length, that is, it damps out at infinite time. One 
example of infinite-length wavelet is 

(1, 1/2, 1/4, 1/8, 1/16, 1/32, 1/64, 1/128, 1/256, •*'). 

3 . Convolution in Discrete-Time Systems 

In order to find a linear system's response or output 
resulting from a specific input signal, we have two methods 
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in general use, One method is to model or characterize 
the linear system by means of linear differential 
(continuous-time systems) or difference (discrete-time 
systems) equations. The output of the system to any input 
signal is then obtained by solving the appropriate equation. 
Another method is based on the superposition property of 
linear systems, using the convolution operation between the 
input signal and the system's impulse response to get the 
output signal. Since we are only interested in treating 
signals in the discrete-time domain, we consider the convo- 
lution operation here only in discrete-time systems. 

The convolution of two wavelets a^ and b is denoted 
by 

c = a * b 



and is defined by the formula 



c. = E a b, 
t m _ A m t-m 
m=0 



To illustrate the use of this formula let us suppose 



a 



= (a Q , a x> a 2 ) = (1, 3, 2) and 



b = (b Q , b 1 , b 2 ) = (2, 1, 4). 



Using the convolution formula, we have 



c a * b ( Cq , c^ , c 2 , Cg , c 4 ) 



where 



c 0 a 0 b 0 1x2 2 



C 1 = a 0 b l + a l b 0 “1 x 1 + 3x 2 =7 



12 



Ix4 + 3xl + 2x2 
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c 2 = a 0 b 2 + a i b l + a 2 b 0 = 



c 3 = a l b 2 + a 2 b l “3x4+2x1=14 



C 4 = a 2 b 2 = 2 x 4 = 8 
Hence 

c = a * b = (2, 7, 11, 14, 8). 

In general case, if a^ and are not wavelet but any two 
digital signals, the formula for convolution operation is 

c = a * b 
where 

OO 

c. = 2 a b , . 

t m t-m 

m=-°° 

4 . The Z-Transform 

In analyzing linear systems, we can use linear 
differential or difference equations. We can also use 
convolution operations. Both are called time-domain tech- 
niques since the independent parameter is time t. But we 
find that if we transform the time signals within our system 
to another form, we can often express relationships more 
simply than we have previously in the time-domain. We call 
this transform approach transform-domain techniques. One 
important consequence of the transform-domain approach of 
linear systems is that the convolution operation in the 
time-domain is converted to a multiplication operation in 
the transform-domain. There are three transforms which are 
important in the study of linear systems. The Laplace 
transform and Fourier transform are used in the continuous-time 



13 



systems. The Z-transform is appropriate for discrete- 
time systems. We will only consider the Z-transform here. 

The Z-transform of finite wavelet is the polynomial in Z 
whose coefficients are the coefficients of the wavelet. For 
example, the Z-transform of the wavelet 

a — (Uq, a.^ , a 2 , . . . , 
is the polynomial 

A(Z) = a^ + a^Z + a 2 Z 2 + , . . + a n Z n . 

The Z-transform of an infinite wavelet is the power 
series in Z whose coefficients are the coefficients of the 
wavelet. For example, the Z-transform of the wavelet 

b = (bQ, b^, b 2 ...) is the power series 

B(Z) = b Q + b 1 Z + b 2 Z 2 + 

The Z-transform of a general two-sided digital signal 
which is of infinite length in both directions will be a 
Laurent series in Z. For example, the Z-transform of the 
digital signal 

c (..., c_ 2 , c_.^ , ^0* ^1* c 2 
is the Laurent series 

C(Z) = . . . + c_ 2 Z“ 2 + c^Z -1 + c Q + c ± Z + c 2 Z 2 + . . . . 

To show convolution operation in time-domain may be 
performed by multiplication operation in transform domain, 
let us redo our example in the last section. We have 

c = a * b 
where 
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a = (a Q , a.^, a 2 > = (1, 3, 2) and 



b = (b Q , b ±> b 2 ) = (2, 1, 4), 



The Z-transforms of a^_ and b^ are 



A(Z ) = 1 + 3Z + 2Z 2 and 
B(Z) = 2 + Z + 4Z 2 . 



Multiplying A(Z) and B(Z), we obtain 



C(Z) = A(Z)B(Z) = 2 + 7Z + 11Z 2 + 14Z 3 + 8Z 4 



After we transform back to time-domain, we get 



c ( Cq , , c 2 , ^ > H > 14, 8) 



which is the same result as in the last section. 

5 . Autocorrelation and Cross-Correlation 

Autocorrelation and cross-correlation are two 
important concepts in studying stationary stochastic processes. 
For our purposes here we need only consider these concepts 
as applied to discrete energy signals. 

The autocorrelation function of any energy signal b^ 
is defined as: 



where the integer x is the time-shift and b^. is the complex- 
conjugate of b^ , If we only consider real-valued discrete 
energy signals, the formula is simply: 



CO 





CO 





When x = 0, we have 
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<f>bb (0) 



E b 
t=-°° 



2 

t 



which is equal to the energy E of the signal, 
If we let s = t + x in the formula, we have 



<P 



bb 



(t) 



oo 


oo 


oo 


E b t + T b t - 


Ebb 
_ s S-T 


= e 


t=-°° 


s = — °° 


s= 



b 



s 




(-T ) . 



Thus the autocorrelation function for a real-valued discrete 
energy signal is an even function of the time-shift t. 

As a generalization of the autocorrelation function, 
the cross-correlation function of two real-valued discrete 
energy signals a, and b. is defined as: 



W 0 = * a t + T b t = E ^S-T 

t=-°° s=-°° 



From the above formula it is easily seen that: 



+ab< T > - 



♦ba ( - T) 



Also we see that the autocorrelation function 4>bb^ T ^ 
a special case of the cross-correlation function 
occurring when the signals a^ and b^_ are the same energy 
signal . 

6 . Deconvolution 

Deconvolution is the inverse operation of convolution. 
Consider the problem that 

c = a * b 



which means the discrete signal c^ is equal to the discrete 
signal a^ convolved with discrete signal b^ . If we are given 
c and a^ and we wish to find b^ , we need to deconvolve a^ 
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from Cj. to get b^, We can do the deconvolution both in 
time-domain and in transform-^domain , 



C t = 



r* 



m=o 



a b^_ 
m t-m 



The first few terms are 



c 


= a b 






o 


o o 






C- 


= a b. 


+ 


a. b 


1 


o 1 




1 o 


Co 


= a b 0 


+ 


a, b_, + a„b 


2 


o 2 




11 2 o 



In the first equation, we can solve b Q , because we know 
both c and a . After we get b , we can use the second 
equation to solve b^. We can continue this iterative process 
until we find the signal b^ , 

In the transform-domain, we know the convolution operation 
becomes a multiplication operation. Hence the deconvolution 
operation becomes a division operation. Thus if c = a * b, 
then after we take the Z-transform, we have 
C(Z) = A(Z)B(Z) , 

If we are given c^ and a^_ , we know C(Z) and A(Z), We can 
find b^ by first solving B(Z) which is equal to C(Z) divided 
by A(Z). Then we can transform B(Z) back to the time-domain 
to get the signal b^ . 

To show the procedures stated above, let us do an 
example. Suppose 
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c t ( Cq , , Cg , Cg ) (12, 10, 14 , 6) 

a t = ( a Q , a x ) = (4, 2). 

We want to find such that 
c = a * b. 



Taking Z-transf orms of and a^ , we obtain 
C(Z) = 12 + 10Z + 14Z 2 + 6Z 3 
A(Z ) = 4 + 2Z 



The Z-transform of b^ 
B(Z) = C(Z)/A(Z) = 



is 

12 + 10Z + 14Z 2 
4 + 2Z 



+ 6Z 3 



3 + Z + 3Z 2 . 



Therefore 



b t = (b 0 ’ b l» b 2 } = (3> 1 > 3) * 

One thing we must point out here is that when we convolve 

two finite-length discrete-time signals, the result is always 

a finite-length discrete-time signal. But when we deconvolve 

a finite-length signal from another finite-length signal, 

the results are not always of finite-length, For example, 

if 



c = a * b 

C t = <V c l» c 2’ C 3 ) = ( ' 12, 10, 14, 8) 
a t = ( a Q , a t ) = (4, 2) 
we will obtain 



b t = <b 0 , b x , b 2> 

= (3, 1, 3, 1/2, - 1/4, 1/8, - 1/16, ...) 
which is of infinite-length. This fact is related to the 
design of least error energy shaping filters, a subject we 
will discuss later on. 
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7. Stationary and Ergodic Stochastic Processes 



A stochastic process is a proces that assigns, 
according to a certain rule, a time function (continuous or 
discrete) f(t, £) to every outcome £ of an experiment. It 
can be viewed as a function of two variables t and c. 
f(t, c) may represent: 

(1) A family of time functions (t and £ variables) 

(2) A single time function (t variable, r, fixed) 

(3) A random variable (t fixed, L, variable) 

(4) A single number (t fixed, r, fixed), 

A stochastic process f(t) (we usually omit its 
dependence on £; ) is stationary in the general sense if its 
expected value is a constant and its autocorrelation depends 
only on t^-tg: 

E{f(t)} = constant ( for any time t) 

E{ f ( t+x )f ( t ) } = R(t) (for any time t). 

A stochastic process f(t) is ergodic if its time averages 
equal its ensemble averages. 

We will assume the earth tilt data as an ergodic 
stationary time series in the design of single-channel 
predictor . 
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B. THEORY OF SINGLE-CHANNEL PREDICTOR 



1 , Method of Approach 

a. Inelegant Method 

By using the FFT (Fast Fourier Transform) we can 
determine the spectral density function of a time series. 
After we invert the spectral density function to the time 
domain we can get a Fourier series representation of the 
time series. We then extrapolate to future time to obtain 
a crude singe-channel prediction. 

b. Elegant Method 

The Wiener-Hopf equation can be used to design 
an optimum filter for a single channel predictor, The 
solution of the Wiener-Hopf equation is facilitated on a 
high-speed digiter computer by using the Levison recursive 
method which can be programmed. We will discuss this method 
of prediction in the following sections. 

2 , Block Diagram of a Linear Time-Invariant Digital 

Filter 

The following is a block diagram of the input/output 
relationship of a linear-system: 




x is the input signal to the linear-system, y is 
the output signal from the linear-system. The relation 
between x and y is: 
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y = x * a 

where a is the impulse response sequence of the linear-system. 

The above equation is the time-domain relation, which 
says the output y is equal to the convolution of the input 
x and impulse response sequence a of the system, 

In transform domain, the input-output relationship 
simply becomes: 



Y(Z ) = X(Z)A(Z) . 

That is, the Z-transform of the output signal y is equal to 
the product of the Z-transform of the input signal x and the 
Z-transform of the impulse response sequence a of the system. 
A system that generates its output in this way is called a 
linear time-invariant digital operator, or a linear time- 
invariant digital filter. The coefficients of the impulse 
response sequence are called the operator coefficients or 
filter coefficients. 

3 . Least-Squares Shaping Filter 

From section 2.2, we know in a linear time-invariant 
system if we are given the input signal x and the impulse 
response sequence a of the system, we can get the output 
signal y by convolving x and a. On the other hand, if we 
are given the input signal x and we want the output signal 
y to be some desired form, how should we choose the impulse 
response sequence or filter coefficients a? Clearly this is 
a deconvolution problem. In transform-domain it becomes 
a division problem. That is, 



A(Z) = 



Y(Z) 

X(Z) 
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where A(Z) is the Z-transform of the desired filter coeffi- 
cients a. But we know from elementary algebra that even Y(Z) 
and X(Z) are both of finite terms, A(Z) may not be always 
of finite terms. In fact in most cases A(Z) will be a 
polynomial of infinitely many terms. This is to mean we 
want a filter with infinite-length. If in the real-world 
situation we can only build a finite-length filter, then we 
can not obtain the exactly desired output signal y. The 
difference between the actual output and the desired output 
is the error. For a specific finite-length filter, the optimum 
one is that one which will yield an output with smallet error. 
Like the regression analysis in statistics, here the smallest 

error means the least squared error. If we denote the desired 

✓\ 

output signal to be y and the actual output signal to be y, 
we want to minimize 



n 

Z 

t=0 



, A x2 

- y t ) 



The optimum filter which can do this job is called a least- 
squares shaping filter since it can generate an output with 
least squared error to the desired shape of output. 

4 . Discrete-Time Wiener-Hopf Equation 

In this section we will discuss how to determine 
the least-squares shaping filter's coefficients. Let the 
input signal be: 

x = (x Q , x ± , • • • , x n ) 

let the desired output signal be 



y = (y 0’ y l’ 



n+m 



) 
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let the filter coefficients be 



f < f 0’ f l’ 



, f ) 
m 



The actual output signal y is the convolution of the input 
signal x and the filter coefficients f, that is: 

y = x * f, 

Now we want to choose f which will minimize 



n+m „ n+m m 9 

I = 2 (y t - y t ) = Z (y - E f X r, 

t=0 z z t=0 x s=0 s r s 



The quantity I is minimized if its partial derivatives 



with respect to each of the coefficients f^, f^, 
equal zero. Therefore we want: 



m 



Si T 


n+m 


m 


d I 

d± r 

3 


o 
w II 

>P 


2(y t - J 0 f s x t-s>(- x t-j> - 0 


(j = 


0, 1 


, 2, • • • , m) 



simplifying, we obtain: 



n+m 



n+m m 



- Z y x . + Z ( Z f x ) x . = 0 
t=0 t t=0 s=0 s t S t " 3 



( j = 0, 1, 2, • • • , m) 



or 



m 


n+m 




n+m 


Z f 


Z 


X, X, 


Z 


s=0 s 


O 

II 


t-s t-J 


o 

II 

•p 


(j = o 


, i, 


2, • • • , m) 





t t-j 
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If we use the notation 



n+m 

r . = E x 



3 ~ s 



x. 



t=0 



t-s t-j 



and 



n+m 



g j 



E y, x, . 
y t t-j 



t=0 



where r j_ s is the autocorrelation (for index j-s) of the 

input signal x and g. is the cross-correlation (for index j) 

0 

of the output signal y and the input signal x, we obtain 



m 



2 f r . 0 = g. (j = 0, 1, 2, m). 



s=0 



s J-S 



After writing out this set of m+1 equations, we have 



i o r o 


+ 


f i r -i + 


V-2 + 


+ f r = g„ 

m -m & 0 


f 0 r l 


+ 


f l r 0 + 


f 2 r -l + 


+ f n r l-m = *1 


Hj 

o 

to 


+ 


f l r l + 


f 2 r 0 + ' ' ' 


+ f m r 2-m = g 2 



f n r + f-r + f 0 r 0 + 
0 m 1 m-1 2 m-2 



+ f m r 0 = V 



or in matrix form 
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— — 




r— - 







r 0 r -l r -2 ’ " ' r m 




f o 




g 0 


r l r 0 r -l ' * ' r l-m 




f l 




g l 


-r r r « • • t 

2 1 r 0 2-m 




f 2 





g 2 


• 




• 




• 


• 




« 




• 


♦ 




t 




• 


r r -i r o • • • 

m m-1 rn-2 0 




f 

m 




g m 


i— — i 


— _ 








The above set of equations is called the normal equations 

and is the discrete-time Wiener-Hopf equation. 

5. A Computer Program for a Single-Channel Predictor 

In a least-squares shaping filter system, if the 

desired output signal y is equal to the input signal x 

shifted t time units ahead, that is 
o ’ 




then the filter becomes a predictor since its objective is 
to predict the value of x^ t Q time units ahead. In this case, 
the discrete-time Wiener-Hopf equation 
m 

2 f r . = g ( j = 0, 1 , 2, • • • , m) 

s=0 s 3 ~ s 3 

simply becomes 

m 

1 r i o = r i +t (3 = 0, 1, 2, •••, m) 

s=0 s J_s 3 o 



because 



n+m 



g i = A y t x t- 



t=0 



t-j 



n+m 

E X t+t 
t=0 X o 



x 



t-j 



r . 

3 



+ t 

o 
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Writing out the equations in matrix form, and by noting 

r = r for a real-value signal, we obtain 
x -t t. » 



r o 


r l r 2 


r m 




~ f o~ 




<-+ 

o 

1 


r l 


r o r i 


r m--l 








Pt o +1 


r 2 


V V • « r 

1 0 


r m-2 




f 2 


— 


r t +2 
o 










• 




* 

% 

% 


r m 


r i r 0 > > , 
m-1 m-2 


r o 




f 

m 




r t +m 
o 



Although we must solve the above m+1 equations to find the 



filter 


coefficients 


f f 

O' 1 ’ 


. . « f 

1 m' 


we can take advantage 


of the 


special form 


of the 


autocorrelation matrix 




r A 


r^ 


r 0 


• • r 






0 


1 


2 


m 






r i 


r o 


r i 


• • r 

m-1 




R = 


r 0 


r., 




• • r 0 




m 


2 


1 


0 

• 


m-2 






r 

m 


4 

3 

1 

h* 


• 

r 

m-2 


•' r o 




in order to 


reduce 


the computational 


work . 



The autocorrelation matrix has the property that all 
the elements on any given diagonal are the same. It is also 
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a symmetric matrix. Hence the entire autocorrelation 

matrix R of order m+1 does contain only m+1 distinct 
m 

elements, namely r^, r^, r^, r m « Let us now see how 

we can take advantage of this structure. 

The method we use is called the Levinson recursive 
method. Suppose that we want to solve the following matrix 
equation : 



— 










— — 




— 


r o 


r l 


r 2 


r 

rn 




f o 




g 0 


r l 


r o 


r l ••• 


3 

1 




f l 




g l 


r 2 


r l 


r o 


3 

1 

DO 




f 2 


= 


g 2 






* 

* 






• 




♦ 


r m 


r m- 


l r m-2 * ‘ * 


r o 




f 

m 




g m 


— 






— 








_ _l 



where the autocorrelation coefficients r^, r^ , r^, ' ‘ » r m 

as well as the cross-correlation coefficients g^, g^ , g 2 , 

• • • , g m represent the known quantities, and the filter 
coefficients f , f ^ , f^, • • • , f m represent the unknown 
quantities. The recursive procedure solves the equation in 
a step-wise manner. The step k=0 is given as an initial 
condition, and then the steps k=l, k=2, •••, k=m are done 

successively in a recursive manner. The desired filter 
coefficients are the ones which result on the completion of 
the final step, namely step k=m. 
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Let us now show how we proceed from step k=n to step 
k=n+l. The completion of step k=n requires that we have 
computed the following quantities: 



a 



nO’ nl' 



i ; a , 
nn ’ n ’ 



n 



and 

^nO’ ^nl’ ^nn’ ^n* 

By definition, these quantities satisfy the matrix equations 



r 0 r l r 2 ’ ' ' r n+l 




a A 

nO 




a 

n 


r l r 0 r l r n 




a i 
nl 




0 


r 2 r l r 0 ’ " r n-l 

♦ 

• 




« 

a nn 


— 


• 

0 


r r r .. r. 

n+1 n n-1 0 




0 




_ 6 n_ 



and 



r 0 r l r 2 * ’ ' r n+l 








g 0 


r l r 0 r l r n 




f nl 




g l 


r 2 r l r 0 * ’ * r n-l 

• 

t 




• 

• 

• 

f 

nn 


— 


g n 


r , - r r - • • • r~ 

n+1 n n-1 0 




0 




Y n 



( 1 ) 



( 2 ) 
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Because of the special structure of the autocorrelation 
matrix, equation (1) may be written in the following 
equivalent form: 



r o 


r l 


r 2 


• • • r 


r l 


r o 


r l 


• • « 


r 2 


r l 


r o 


• • • p 






t 




r n+l 


r 

n 


r n-l 


• • • p 



n+1 

n 

n-1 



0 



0 



a 



nn 



a 



nl 



nO 



n 



0 



0 



a 



n 



(3) 



Let C be a constant which we will determine in order to get 
a desired result. The desired result is after we multiply 
equation (3) by C n and add the result to equation (1), we 
want it to be identical to a new equation. That is, we want 
equation (4) and (5) to be identical. 



r o 


r l 


r 2 


« • • yi 

r n+l 




a n0 


+ 


0 




a 

n 


+ 


C 3 
nn 


r l 


r o 


r l 


... r 

n 




a nl 


+ 


c a 
n nn 






0 




r 2 


r l 


r o 


r n-l 






* 

• 




= 




* 














a nn 


+ 


c a ^ 
n nl 






0 




r n+l 


r 

n 


r n-l 


r o 




0 


+ 


c a n 
n nO 




b „ 


+ 


c a 
n n 



(4) 
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r 0 r l r 2 • ' • r n+l 




a n+l ,0 




a , - 
n+1 


r l r 0 r l * ” r n 




a n+l , 1 




0 


Y Y Y * • . r* 

2 1 0 n-1 

* 




* 

• 




t 

• 


* 




« 

a n+l , n 




t 

0 


r n+l r n r n-l ‘ * ' r 0 
L 




a n+l , n+1 




0 



In order for equations (4) and (5) to be identical, we 
require that 



a 



n 



c 3 
n n 



a n+l 

0 



0 



3 + 

n 



c a 
n n 



0 

0 



( 6 ) 



and 



p 1 
p 
O 

+ 

o 

_ . J 




a n+l , 0 


a - + c a 
nl n nn 

* 




a n+l , 1 

t 


• 

* 

a + c a - 

nn n nl 




% 

* 

a n+l , n 


0 + c a _ 

n nO 




a n+l , n+1 



30 



From equation (6), we can determine the value of c n by 
6 n + °n a n = 0 



or 



n 



n a 



n 



After has been determined, we can compute all the 



following quantities : 



a n+l ’ a n+l , 0 ’ a n+l,l’ a n+l,2’ ”*» a n+l,n+l* 

Next, let 

^n+1 - a n+l , O r n+l , 0 r n+2 + a n+l,l r n+l + + a n+l,n+l r l 



then we get the counterpart of equation (1) of step k=n in 
step k=n+l as follows: 



r 0 r l r 2 ' " r n+2 




a n+l , 0 




a . 

n+1 


r l r 0 r l r n+l 




3- , i i 

n+1 , 1 




0 


r „ r- r~ • • • r 










2 10 n 




t 




« 


t t • 




• 

a n+l , n+1 




0 


r n+2 r n+l r n r 0 




0 




g n+l 



Now, rewrite equation (5) as follows; 
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r 0 r l r 2 ' ' ’ r n+l 

r l r 0 r l * " r n 

r 2 r l r 0 ' ’ ‘ r n-l 




a n+l , n+1 
a 

n+1 , n 
• 

• 

• 




r 

o o - • • 

1 


« 

r r r - • • • 

n+1 n n-1 0 




a n+l , 1 
a n + l , 0 




0 

a . -i 
n+1 



Let us multiply equation (9) by a constant d n , as yet under- 
undetermined, and then add the result to equation (2). We 
obtain 



r 0 r l r 2 • • • r n+l 




f n0 +d n a n+l,n+l 




V° 


r, r n r- • • • r 




f .. +d a , - 




g-i + o 


10 1 n 




nl n n+l,n 




b l 


r 2 r l r 0 ’ ’ ' r n-l 




• 

• 


_ 


• 


• • • 




f +d a - 

nn n n+i , 1 




V° 


r n+l r n r n-l * * * r 0 




0 +d a - _ 
n n+1,0 




Wn+1 
L_ _J 



We want equation (10) to be identical to the following 
new equation: 



r 0 r l r 2 * * * r n+l 




f n+l,0 




i 

o 

fcJD 

1 


r l r 0 r l * * * r n 




f n+l, 1 




g i 


r 2 r l r 0 * ’ * r n-l 

% 




« 

« 

f L 1 
n+1 , n 


— 




r ,i r r - • • • r A 

n+1 n n-1 0 




f 

n+1 , n +1 




_ gn+ L 



(ID 
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This means we require 



g o 




g 0 


g l 




g l 


« 






* 


= 


• 


* 




t 


g n 




g n 


’ Wn+l 




g n+l 


— 




l 



and 



( 12 ) 











f n +d a , 
nO n n+1, n+1 




f n+1,0 


f .. +d a ... 
nl n n+l,n 

• 




f n+l,l 

• 


f +d a , 
nn n n+1 , 1 




t 

• 

f 

n+1 , n 


d a ... _ 
n n+1,0 




, f 

n+1 , n+1 



From equation (12), we can determine d R by the formula 
d . 



After d R is substituted in equation (13), we can compute the 
following quantities: 

f n+1 , 0 ' f n+1 , 1 ’ f n+1 , 2 ’ ' * * ’ f n+1 , n+1 * 
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Next, let 



Y n+1 ^n+1 , 0 r n+2 + ^n+l , l r n+l + ' ' * + ^n+l , n+l r l 

We obtain the counterpart of equation (2) of step k=n in 
step k=n+l as follows: 



r 0 r l r 2 r n+2 




f n+l,0 




“T1 

(JQ 

O 

1 


r l r 0 r i r n+l 




^n+l , 1 




g l 












2 10 n 








] 


• 




• 




f 


• 




^n+l , n+1 




g n+l 


r . o r , - r • « • r_ 
n+2 n+1 n 0 




0 




Y n+1 



We have completed the recursive procedure from step k=n to 
step k=n+l . To start at step k=0, we must determine 

a 00 ; a 0’ e o 

and 



f 00 ; Y 0 

which satisfy the matrix equations: 



r o 


r l 




a oo 




a o 


r l 


r o 




0 







and 



( 15 ) 
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( 16 ) 



r o 


r l 




f oo 




g 0 


_ r l 


r o_ 




0 




Y o 



In equation (15), we set a Q0 = 1 and solve 
we obtain 



a 0 r 0’ 



From equation 
Sr 



00 



’0 

'o 




(16), 



Y 



0 



1' 

we obtain 



= v f 

1 00 ' 



Starting from step k=0 as above, we can proceed recursively 
to steps k=l , 2, 3, •••, m. At the final step k=m, the 

values obtained for the filter coefficients, namely 



f f f . • • f 
mO ’ ml’ m2 * ’ 1 mm 

will be the required solution. A computer program based on 
the above recursive method is included in Appendix A. 
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Ill, CHARACTERISTICS OF TIDAL DATA 



A. INTRODUCTION 

By Newton's law of universal gravitation, each body in 
the universe attracts another body with a force directly 
proportional to the product of their masses and inversely 
proportional to the square of the distance between them. 

We may express this fact in a formula: 

Gm,m n 



where is the gravitational force, and m 2 are the 

masses of the two bodies in the universe and r is the distance 

between them. G is called the universal constant of grav- 

-11 2 2 

itation and is approximately equal to 6.67 x 10 n*m /kg . 

Celestial bodies such as the earth, the sun, and the 
moon are, of course, obeying the law of gravitation. Because 
they are not point masses and have large dimensions, different 
points on the surface of the body will be affected by differ- 
ent forces. Tidal forces are generated due to this effect. 

In this chapter we will give a general description of the 
characteristics of the earth tidal data. 

B, GRAVITATIONAL POTENTIAL AND TIDAL POTENTIAL 
1 . Basic Formulas From Mathematics 

Since in the derivation of tidal potential we will 
use some formulas from applied mathematics and spherical 
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trigonometry, we give a brief discussion of them in this 
section . 

a, Legendre Polynomials 

The Legendre Polynomials are defined by: 

p _ ( 2n-l ) ( 2n-3 ) . . . 1 f„n n(n-l) „ n-2 

n tx; n! j_ x " 212^17 x 

. n(n^-l)(n-2)(n-3) n-4 

2.4(2n-l)(2n-3) 

We observe that P n (x) is a polynomial of degree n. The first 
few Legendre polynomials are as follows: 




P Q (x) = 1 P 3 (x) = l/2(5x 3 -3x ) 

P x (x) = x P 4 C X ) = l/8(35x 4 -30x 2 +3) 

P 2 (x) = l/2( 3x 2 -1 ) P 5 (x) = l/8(63x 5 -70x 3 + 15x). 



The Legendre polynomials can also be expressed by Rodrique's 
formula : 



,n 



V x > * 7H- 

2 n ! 



dx 



n 



, 2 - .n 

(x - 1) . 



After we have obtained Pq( x ) and P^(x) by the formulas 
given above, we can also use the following recurrence 
formula to obtain the higher degree Legendre Polynomials: 



P n + l (x) = 



2n-t-l 

n+1 



xP (x) - 
n v 7 



n 

n+1 



P n-! (X >’ 



The function is called the generating function 

hj l-2xt+t 2 

for Legendre polynomials. By using the binomial theorem 
we can prove the following result: 
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1 



oo 



= Z P (x)t n . 
n=0 n 

b. The Cosine Law for Spherical Triangles 

In a spherical triangle, if we denote the three 
sides to be a, b, c and the three angles to be A, B, C; the 
following relations hold. 

cos a = cos b cos c + sin b sin c cos A, 



cos A = - cos B cos C + sin B sin C cos a. 



We call the above formulas the cosine law for spherical 
triangles . 

2 . Derivation of Gravitational Potential and Tidal 
Potential . 

a. Gravitational Force, Gravitational Field 
Intensity, and Gravitational Potential 

According to Newton's law of universal gravita- 
tion, a point mass m^ at a distance r from another body of 
mass M will experience a force 



F = -G 
g 



m 0 M 



where i is a unit vector in the direction of r. F is 
r g 

called the gravitational force. The minus sign denotes it 
is an attractive force. 

We define the gravitational field intensity at a 
point to be the gravitational force experienced by a unit 
mass at that point. If the unit mass is at a distance r 
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from the body with mass M, we will have 



i- - -g -M r 

g 2 r 



where I denotes the gravitational field intensity. 

o 

The gravitational potential V at a point which 

g 

lies at a distance r from a body with mass M is defined as 
follows : 



oo oo 

- f r . d? = f -G — g i- 

g J r g r 2 ] 



dr 



(X) 

- J -G-S 



9 dr = -G — . 
2 r 



Another way to show the relation between V and I is 

g g 



I = -VV = -grad V 
g g g 



Since V is a scalar while F and I are vectors, V is 
g g g g 

-)■ -)■ 

much easier to handle than F and I , In deriving the 

g g 

theoretical tide on earth induced by the moon and the sun, 

we will use tidal potential instead of tidal force. 

b. Gravitational Potential and tidal potential at 
a Point on the Surface of Earth 
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B 



Figure 1 



In Figure 1, let E be the earth, P is any point 



on the surface of E. The gravitational potential at P created 
by a celestial body of mass M is 




where R is the distance PB. Notice we are assuming that the 
whole mass of the celestial body is concentrated at a point 
B. By the cosine law for plane triangles, we have: 

D 2 2 2 0 _ 

R = a + r - 2arcos 

Therefore we can rewrite the gravitational potential V as 



follows : 



V = -G 



M 





- 1/2 
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After we expand the above equation, we get a more compact 
form : 



,, 00 n 

v = -G - Z P (cos O(-) 
r ~ n v 
n=0 

where P n (cos £ ) are the Legendre polynomials. The first 
term of V is 




which is a constant and has no contribution to the gravitational 
force. Therefore it has no physical significance and we will 
disregard it, The second term is 

V 1 = -G - (cos C) — = -G — 2 a cos £ . 

r 

To evaluate the contribution of V 1 to the gravitational 
force, let us choose a cartesian frame of reference centered 
at 0, the x axis running in the direction OB and the y axis 
running in the direction 00'. With this choice of coordinate 
system we have 

a cos £ = x , a sin ( = y , 

and the contribution of to the gravitational field intensity 
is 

!+ = _ VVl = -V( -G - 2 x) = G -f i* 

r r 

where i** is a unit vector in the direction of x axis. This 
x 

shows that 1^ is a uniform gravitational field. It has a 

M 

constant magnitude G and a constant direction which is 

r 

parallel to the x axis. Hence will not create tidal effect. 
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The third term of V 



is 



V 



2 



-G “ • 1/2 

-G ^3 • 1/2 
r 

-G ^3 • 1/2 
r 

-G ^3 • 1/2 



9 2 
(3 cos^C - 1) (f) 

, 0 2 2 2 
(3a cos £ - a ) 

/o 2 2 2 N 

( 3x - x - y ) 

/o 2 2, 

( 2 x - y ) 



in our choice of coordinate system. The gravitational field 
intensity created by is 



->■ 





-V 






) 



= G 



M 



0 ( 2 xi -yi) 
3 v x J y 



where i and i are unit vectors in the direction of x axis and 
x y 

y axis. We observe that is different at different locations 
because it depends upon x and y. This non-uniform gravita- 
tional field intensity will create a tidal force at all points 
on the earth. Hence i s responsible for the tidal effect. 

We can also show that after Vg, every term of V is responsible 
for the tidal effect. Therefore we call 

M 1=0 o n 

Vrp = -G — E P (cos £) (-) 

T r n v r' 

n — 2 



the tidal potential. 

3 . Tidal Potential on the Surface of Earth Induced by 
the Moon and the Sun 

Although theoretically every celestial body in the 
universe creates a tidal potential on our earth, they are so 
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small that in practice we can neglect all of them except 

those created by the moon and the sun. If we denote the tidal 

potential at a point P on the earth induced by the moon be 

V , by the sun be V , we have 
m’ J s' 



m 



= -G — 



M 

_r 

r 



n 



m 



Z P (cos ? ) (" ) 
n=2 m 



V 



s 




Z P (cos £ ) ( — ) 
n=2 s 



where 

22 

M = the mass of the moon - 7.38 x 10 J kg, 
m ’ 

30 

M = the mass of the sun - 1.991 x 10 kg, 
s 

r = the distance from the center of the earth to the 
center of the moon, 

r = the distance from the center of the earth to the 
center of the sun, 

£ = the zenith angle of the moon, 

m 

£ = the zenith angle of the sun, and 

s 

a = the radius of the earth - 6365 km. 



If we use the mean value of r and the mean value of r to 

m s 

3. cl 

calculate — and — , we obtain 
r r ' 



a jt 6365 
r m 3*84xl0 5 



1 • 66 x 10 




6365 

1 • 49xl0 8 



4*33 x 10^ 5 
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cL cL 

Since — - and — are both very small, we can express V and 
r r m 

m s 

V approximately as follows: 
s 

M 2 

V ± _g — P 0 (cos £ ) ( — ) 
m r 2 v ^m' r ' 

m m 



= -G 



M a 
m 



m 



• 1/2 ( 3cos K -1), 
' m 



M 2 

V = -G — P 0 (cos £ ) ( — ) 
s r 2 v s' v r ' 

s s 



= -G 



M a 
s 



1/2 ( 3cos C s -1). 



For the same value of C and £ , the ratio is 

m s ’ 




This indicates that the tidal potential due to the moon 
is roughly twice that due to the sun. Therefore the moon 
should be the preponderant source of the tidal effect on 
our earth, as is indeed observed. 
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C. THE PERIODICITIES IN THE TIDAL POTENTIAL INDUCED BY THE 
MOON AND THE SUN 



CELF.ST 



POLE 




ZENITH 




Figure 2 
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In Figure 2, let P be a point on the earth with latitude 

L and longitude A , Q be the geographical point of the 

celestial body M which has declination d and local hour angle 

t at the instant shown in the figure. If the celestial body 

is the moon, then by the result of the last section we can 

write the tidal potential at P induced by it as: 

„ 2 

M a 9 

V - -G • 1/2 (3 cos" a-1). 

m o' 

r 

m 

Now we will derive the periodicities present in the tidal 

potential. By using the cosine law for spherical triangles, 

we can rewrite V as 

m 



M a* 

V - -G - 51 o- 
m 3 

r m 



= -G 



M a" 
m 



m 



1/2 J^3( cos Ocos Z + sin 6sin Z cos t) 2 ~l] 

1/2 |j}(sin L sin d + cos L cos d cos t 2 -lj 



After we expand the above equation and algebraically 
simplify, we obtain 



V = -G 



m 



M a‘ 
m 



m 



3/4 [^3(sin 2 L - 1/3) (sin 2 d - 1/3) 



2 2 

+ sin 2L sin 2d cos t + cos~L cos~d cos 2t 



[3(si: 



‘] 



K m | 3(sin 2 L-l/3) (sin 2 d-l/3] 



9 2 

+ k (sin2L sin 2d cos t) + K (cos~L cos d cos 2t) 
m m 



where 
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is a constant. 



K = -3/4 G 
m ' 



M a 
m 



m 



L is the latitude of P and is also a constant. d and t are 

the declination and local hour angle of the moon, respectively 

They are changing from time to time. t will go through a 

complete cycle of 360° in one lunar day which is about 24 

hours and 50 minutes on the average. d will go through a 

complete cycle in about one month. Since in one lunar day 

d will change very little, we will consider the term 

K(sin 2L sin 2d cos t) as the tidal potential component with 

2 2 

frequency 1 cycle/lunar day and the term K(cos L cos d cos2t) 
with frequency 2 cycle/lunar day. We also call them diurnal 



and semi-diurnal component. The first term K 



3(sin L-l/3) 



2 ~l 2 

(sin d-1/3) contains sin d which is an even function of d. 
Since d goes from zero to positive maximum and then back to 



zero and then goes to negative in a symmetric way 
2 



, K jl 



3(sin L 



-1/3) (sin~d-l/3 )j will be a component with frequency 1 cycle 
per half month. We call it the fortnightly component. 

By the same procedure we can express the tidal potential 
at the point P induced by the sun as 



V = K 
s s 



[3(sin 2 L-l/3)(sin 2 d-l/3)J 



2 2 

+ K (sin2L sin 2d cos t) + K (cos L cos d cos2t) 
s s 



where K 

s 



M a 
s 



-3/4 G -J 
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d and t of the sun go through a complete cycle with periods 

one year and one solar day, respectively. Hence the three 

terms of V have frequencies 1 cycle/half year, 1 cycle/solar 
s 

day, and 2 cycle/solar day. 

If we do not neglect the higher order terms in the tidal 
potential, and also consider the quite complicated relative 
motion of the moon and the sun with respect to our earth 
in detail, the spectrum of the tidal potential should be 
composed of an infinite number of terms. But for practical 
purposes it is quite enough to consider only the most predom- 
inant nine terms, Table 1 is a summary of these terms. 



NAME PERIOD (hr) 



FREQUENCY 
( degree/hr ) 



FREQUENCY (rad. /hr) 



01 



K1 



0 . 23934464E 02 0.15041068E 02 0.26251602E 00 

0 . 25819336E 02 0.13943035E 02 0.24335176E 00 



Ql 

M2 



PI 



M4 



K2 



N2 



S2 



0 . 24065887E 02 0.14958931E 02 0.26108247E 00 
0 . 26868347E 02 0.13398661E 02 0.23385066E 00 
0 . 12420603E 02 0.28984100E 02 0.50586778E 00 
0 . 12000000E 02 0.30000000E 02 0.52359855E 00 
0 . 12658348E 02 0.28439728E 02 0.49636668E 00 
0 . 11967240E 02 0.30082123E 02 0.52503186E 00 
0 . 62103014E 01 0.57968201E 02 0.10117350E 01 



TABLE 1. Summary of Nine Principal Tidal Lines. 
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D. MISCELLANEOUS NOISE IN THE REAL TIDAL DATA 



Meteorological activities such as the rainfall, diurnal 
and seasonal temperature variations and pressure changes, 
introduce noise in the tidal tilt data. Other factors such 
as the ocean loading effect, daily and weekly traffic patterns 
also introduce some noise. However, to a lesser degree, 
these effects can be considered stationary and predictable 
over short periods of time. The effects due to local earth- 
quake activity or underground nuclear tests would introduce 
excessive noise into the real tidal tilt data. These noises 
are random in nature compared to other noise. We will 
remove these noises from the real tidal data when we use 
them as the input data of the single-channel predictor. 



49 



IV. APPLICATIONS OF SINGLE-CHANNEL PREDICTOR 



A. DISCUSSION OF THE PARAMETER’S VALUE USED IN THE COMPUTER 
PROGRAM 

Before using the "PREDICT" computer program, the first 
step is to convert a stationary continuous time signal to 
a discrete one by the method of sampling. The sampling 
period T should be chosen according to the following formula: 



where f is the highest frequency component of the time 
signal. After we have obtained the discrete time signal, the 
next problem we will consider is how to choose the length 
of the input data; how far ahead should we predict and how 
long the filter length should be. Normally we should first 
do a spectral analyssi of the signal; find out the highest 
and the lowest frequency components. Then we should choose 
the length of the input signal data greater than the longest 
period of the spectral components. For the predicted 
distance, it is a good rule to choose it less than one- 
twelfth of the longest period of the spectral components and 
at the same time close to a small multiple of the shortest 
period of the spectral components. For the earth tidal tilt 
data, we have shown in chapter III that the most important 
spectral components are semi-diurnal, diurnal and fortnightly 
Hence we choose the sampling period as one hour which is 
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considerably less than l/2f m (about 6 hours). Also, we 
choose the length of the input data to be 720 hours which 
is a little more than twice the period of the fortnightly 
components. The predicted distance we use is 24 hours which 
is less than one-twelfth of a fortnight and at the same time 
almost equal to two periods of the semi-diurnal components 
and one period of the diurnal components. 

Another parameter used in the "PREDICT" program is the 
length of the filter coefficients. It is the most important 
parameter which affects the computer CPU time used to do all 
the computations of the program. We will discuss it in 
detail in the next section. 

B. AN INVESTIGATION OF THE FILTER LENGTH USED IN THE 

"PREDICT" PROGRAM 

For a signal of pure sine wave form, it is easy to show 
that a filter length of two is enough to predict the future 
data. From the basic trigonometry, we have 



sin 


(t+T) = sin t 


cos T + 


cos t 


sin 


T 


(1) 


sin 


(t-T) = sin t 


cos T - 


cos t 


sin 


T 


(2) 



By adding equation (1) and (2) together, we obtain 

sin (t+T) + sin (t-T) = 2 cos T sin t (3) 

or 

sin (t+T) = 2 cos T sin t + (-1) sin (t-T) (4) 
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Equation (4) means we can use a two-length filter 



coefficients 

(2 cos T, - 1) 



to predict one unit sampling period T ahead if the signal is 
a pure sine wave. More specifically, let us suppose there 
is a sine wave signal with maximum amplitude 1 and sampling 
at a uniform time interval T which corresponds to 30 degrees. 
The discrete time signal will be 



(• 



0 . 






1 ^3 

2 2 ’ ’ 



£ 

2’ 







- 1 , 



The two-length filter coefficients in this case are, 



(2 cos 30°, -1) = ( -1). 



If we convolve the signal with the filter coefficients, we 
obtain a new signal 



(• 



0 , 




i il i 
2 ’ 2 ’ X ‘ 2 ’ 

t 






- 1 , 



which is just the original signal shifted one sampling period 
ahead. If we want to predict two sampling period ahead, we 
can obtain the required filter coefficients in a recursive 
way as follows: 
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sin ( t+2T ) = 2 cos T sin (t+T) + (-1) sin t 

= 2 cos T [2 cos T sin t + (-1) sin (t-T)] 

+ (-1) sin t 

= (4 cos 2 T-l ) sin t + (-2 cos T) sin (t-T) 

which in our example will be, 

(4 cos 2 30°-l, -2 cos 30°) = (2, - ^3) . 

If we convolve the original signal with (2, - ^3) , we will 
obtain a signal, 




i 



which is the original signal shifted 2 sampling periods 
ahead. 

When a signal's spectrum contains more than one 
frequency, the sampling period T will correspond to different 
angle measurements for different frequency components. If 
we want to predict such signals one sampling period T ahead, 
the filter length in general will not be of finite length. 

With a finite length filter we can only approximately predict 
the time shifted signal. 

The earth tidal tilt data contains nine principal' spectral 
components and various sources of noises. To predict the 
earth tidal tilt data by using the "PREDICT" program, the 
best way to choose the filter length is by doing a statis- 
tical analysis for some samples and then make a decision. 
Appendix B contains a summary of test results of RMS error 
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VS. filter length used in the "PREDICT" program. Seven 
sample runs have been made. The conclusions are as follows: 

(1) The RMS error of the last 24 hours is much smaller 
than the RMS error of the whole 720 hours. 

(2) The filter length at 720, which is the maximum length 
for a 720 input data, causes larger prediction error 
than other filter lengths. 

(3) The filter length between 60 and 660 does not have too 
much effect on the predicted error, with only a little 
improvement in the central region. 

Based on the above results, we think that normally a filter 
length of 120 is enough for the prediction. But if we are 
quite concerned about the accuracy of the prediction, we may 
use a filter length near 360. 

If the length of input data is not 720, the rule of 
choosing filter length will be 0.2 to 0.5 of the input data 
length. The length of the input data should be greater 
than the period of the lowest frequency component of the 
signal . 

C. USING THE SINGLE-CHANNEL PREDICTOR AS AN INTERPOLATOR 
If for some reason there is a missing data segment in 
an otherwise continuous record, we can apply the "PREDICT" 
program to predict the lost data. Since theoretically the 
program can be used to predict in both directions, it is a 
good rule to predict the lost data from both the forward 
direction and the backward direction and then take the average 
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value as the predicted data value. It is easy to show that 
the expected mean squared error of the average value is less 
than or at most equal to the expected mean squared error 
caused by one direction prediction. Suppose at an arbitrary 
point in time the prediction error from the forward direction 
is E^, from the backward direction is , the averaged error 
is (E^ + E 2 )/2, we have, 



1 

2 





V E 2 



Ei 2 + E 2 2 -2 Ei E 2 



- ? < E 1 - e 2> 2 £ °- 



or 



I E 1 2 + I e 2 2 £ 



IE 1 + E 2 



\2 



The above inequality means whenever we apply the single- 
channel predictor as an interpolator, we should always predict 
from both directions and use the averaged prediction value. 
Appendix C contains the results of six test runs on the real 
tidal tilt data to illustrate the application of singe- 
channel predictor. 
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APPENDIX A 



COMPUTER PROGRAM 



'PREDICT' 
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APPENDIX B 



RESULTS OF SEVEN TEST RUNS OF PREDICTION 
UNDER DIFFERENT FILTER LENGTHS 



DATA USED 



1 , 

2 . 



3. 



Instrument Number 07 
Instrument Location 

Berkeley California 



City 

Longitude : 
Latitude : 
Date of Data 
Year 1970 

Year 1971 



State 

122°14 • . 1 W 
37°52 ' ,6 N 
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TEST RUN NO, 1 



FILTER 




LENGTH 


MONTH 


60 


0,405 


120 


0.395 


180 


0.386 


240 


0,380 


300 


0.379 


360 


0.370 


420 


0,367 


480 


0.367 


540 


0 , 365 


600 


0.357 


660 


0.354 


720 


0.349 



RMS ERROR 

LAST DAY 

0,139 

0.132 

0.106 

0,092 

0,089 

0,097 

0.114 

0.117 

0.121 

0.151 

0,143 

0.159 
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TEST RUN NO, 2 



FILTER 


RMS 


ERROR 


LENGTH 


MONTH 


LAST DAY 


60 


0,289 


0.229 


120 


0.281 


0.238 


180 


0,276 


0.255 


240 


0.274 


0.258 


300 


0.273 


0.257 


360 


0.273 


0.264 


420 


0.272 


0.261 


480 


0.271 


0.272 


540 


0.270 


0.275 


600 


0.270 


0.283 


660 


0.270 


0.292 


720 


0.270 


0.302 
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TEST RUN NO, 3 



FILTER 




LENGTH 


MONTH 


60 


0,304 


120 


0.301 


180 


0.299 


240 


0.290 


300 


0.284 


360 


0.282 


420 


0.279 


480 


0.278 


540 


0.278 


600 


0.277 


660 


0.277 


720 


0,277 



RMS ERROR 

LAST DAY 

0.144 

0.135 

0.135 

0.109 

0.131 

0.142 

0.133 

0.127 

0.128 

0.130 

0.131 

0.134 
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TEST RUN NO. 4 



FILTER 


RMS 


ERROR 


LENGTH 


MONTH 


LAST DAY 


60 


0.393 


0.465 


120 


0,407 


0.474 


180 


0.412 


0.472 


240 


0,413 


0.449 


300 


0.413 


0,458 


360 


0.414 


0,463 


420 


0.417 


0.464 


480 


0.422 


0.454 


540 


0.421 


0.464 


600 


0.419 


0.456 


660 


0.422 


0.444 


720 


0.428 


0.615 
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60 

120 

180 

240 

300 

360 

420 

480 

540 

600 

660 

720 



TEST RUN NO, 5 



RMS ERROR 



MONTH 


LAST DAY 


0.194 


0,118 


0.191 


0.113 


0.197 


0.117 


0.200 


0,118 


0.203 


0,119 


0.203 


0.119 


0.205 


0.122 


0.203 


0.120 


0.203 


0.121 


0 . 204 


0.122 


0.205 


0.120 


0.207 


0.144 
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60 

120 

180 

240 

300 

360 

420 

480 

540 

600 

660 

720 



TEST RUN NO, 6 



RMS ERROR 



MONTH 


LAST DAY 


0.409 


0,129 


0.414 


0.136 


0.411 


0.139 


0.410 


0.134 


0.410 


0.126 


0.410 


0.142 


0.407 


0.135 


0.408 


0.131 


0.378 


0.175 


0.375 


0.190 


0.376 


0.186 


0.375 


0.271 
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TEST RUN NO. 7 



FILTER 




LENGTH 


MONTH 


60 


0.368 


120 


0.357 


180 


0.355 


240 


0.359 


300 


0.350 


360 


0.359 


420 


0.367 


480 


0.367 


540 


0.369 


600 


0.369 


660 


0.371 


720 


0.378 



RMS ERROR 

LAST DAY 

0.174 
0.171 
0.153 
0. 168 
0.178 
0.171 
0.165 
0.164 
0.177 
0.194 
0.191 
0.425 
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APPENDIX C 



RESULTS OF SIX TEST RUNS OF DATA INTERPOLATION 
BY APPLYING SINGLE-CHANNEL PREDICTOR 



DATA USED 

1, Instrument Number 07 



Instrument 


Location 




Berkeley 


California 


City 


State 


Longitude : 


122°14 ' 


.1 W 


Latitude : 


37°52 ' 


.6 N 


Date of Data 




April 22, 


1971 




April 30, 


1971 




May 8, 


1971 




May 15, 


1971 




May 24 , 


1971 




May 30 , 


1971 
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April 22, 1971 



ACTUAL 

DATA 

VALUE 

34.950 

34,990 

35.100 

35.260 

35.400 

35.440 

35.380 
35.280 
35.100 
34.900 
34.700 
34.550 
34.450 
34.420 
34.470 
34.580 
34.730 
34.830 
34.850 
34.820 
34.720 
34.580 
34.450 

34.380 



FORWARD 



PREDICTION 



VALUE 


ERROR 


34.996 


-0.046 


35.124 


-0.134 


35.304 


-0.204 


35.450 


-0.190 


35.522 


-0.122 


35.540 


-0.100 


35.502 


-0.122 


35.444 


-0.164 


35.364 


-0.264 


35.250 


-0.350 


35.159 


-0.459 


35.102 


-0.552 


35.100 


-0.650 


35.158 


-0.738 


35.247 


-0.777 


35.340 


-0. 760 


35.418 


-0.688 


35.463 


-0.633 


35.404 


-0.554 


35.277 


-0.457 


35.141 


-0.421 


34.990 


-0.410 


34.894 


-0.444 


34.832 


-0.452 



BACKWARD 

PREDICTION 



VALUE 


ERROR 


34.540 


0.410 


34.631 


0.359 


34.806 


0.294 


34.993 


0.267 


35.153 


0.247 


35.257 


0.183 


35.300 


0.080 


35,255 


0.025 


35.146 


-0.046 


34.994 


-0.094 


34.821 


-0.121 


34.674 


-0.124 


34.580 


-0.130 


34.565 


-0.145 


34.635 


-0.165 


34.769 


-0.189 


34.916 


-0.186 


35.056 


-0.226 


35.163 


-0.313 


35.221 


-0.401 


35.232 


-0.512 


35.190 


-0.610 


35.146 


-0.696 


35.137 


-0.757 



AVERAGE 

PREDICTION 



VALUE 


ERROR 


34.768 


0.182 


34.877 


0.113 


35.055 


0.045 


35.221 


0.038 


35.337 


0.063 


35.398 


0.042 


35.401 


-0.021 


35.349 


-0.069 


35.255 


-0.155 


35.122 


-0.222 


34.990 


-0.290 


34.888 


-0.338 


34.840 


-0.390 


34.861 


-0.441 


34.941 


-0.471 


35.054 


-0.474 


35.167 


-0.437 


35.259 


-0.429 


35.283 


-0.433 


35.249 


-0.429 


35.186 


-0.466 


35.090 


-0.510 


35.020 


-0.570 


34.984 


-0.604 



RMS ERROR FOR FORWARD PREDICTION = 0.464 
RMS ERROR FOR BACKWARD PREDICTION = 0.338 
RMS ERROR FOR AVERAGE PREDICTION = 0.356 
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April 30, 1971 



ACTUAL 

DATA 

VALUE 

33.360 

33.430 

33.500 

33.510 

33.470 

33.460 

33.450 

33.480 

33.530 

33.590 

33.660 

33.690 

33.680 

33.590 

33.420 

33.250 

33.110 

32.990 

32.930 

32.920 

32.920 

32.960 

33.070 

33.170 



FORWARD 

PREDICTION 



VALUE 


ERROR 


33.920 


-0.560 


33.900 


-0.470 


33.890 


-0.390 


33.890 


-0.380 


33.890 


-0.420 


33.930 


-0.470 


34.010 


-0.560 


34.130 


-0.650 


34.210 


-0.680 


34.240 


-0.650 


34.180 


-0.520 


34.080 


-0.390 


33.960 


-0.280 


33.770 


-0.180 


33.540 


-0.120 


33.350 


-0.100 


33.240 


-0.130 


33.180 


-0.190 


33.170 


-0.240 


33.220 


-0.300 


33.330 


-0.410 


33.450 


-0.490 


33.590 


-0.520 


33.680 


-0.510 



BACKWARD 



PREDICTION 



VALUE 


ERROR 


32.860 


0,500 


33.000 


0.430 


33,100 


0,400 


33.120 


0.390 


33.120 


0,350 


33.120 


0.340 


33.080 


0.370 


33.040 


0.440 


33.010 


0.520 


33.000 


0.590 


33.020 


0.640 


33,030 


0.660 


33.070 


0.610 


33.110 


0.480 


33.130 


0.290 


33.110 


0.140 


33.010 


0.100 


32.870 


0.120 


32.760 


0.170 


32,650 


0.270 


32.610 


0.310 


32.590 


0.370 


32.620 


0.450 


32.710 


0.460 



AVERAGE 



PREDICTION 



VALUE 


ERROR 


33.390 


-0.030 


33.450 


-0.020 


33.495 


0.005 


33.505 


0.005 


33.505 


-0.035 


33.525 


-0.065 


33.545 


-0.095 


33.585 


-0.105 


33.610 


-0.080 


33.620 


-0.030 


33.600 


0.060 


33.555 


0.135 


33.515 


0.165 


33.440 


0.150 


33.335 


0.085 


33.230 


0,020 


33.125 


-0.015 


33.025 


-0.035 


32.965 


-0.035 


32.935 


-0.015 


32.970 


-0.050 


33.020 


-0.060 


33.105 


-0.035 


33.195 


-0.025 



RMS ERROR FOR FORWARD PREDICTION = 0.435 
RMS ERROR FOR BACKWARD PREDICTION = 0.421 
RMS ERROR FOR AVERAGE PREDICTION = 0.072 
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May 8, 1971 



FORWARD 
PREDICTION 
VALUE ERROR 

-0.083 



BACKWARD 
PREDICTION 
VALUE ERROR 



AVERAGE 
PREDICTION 
VALUE ERROR 



ACTUAL 

DATA 

VALUE 

32.240 

32.270 

32.340 

32.450 

32,580 

32.660 

32.680 

32.650 

32.550 

32.410 

32.190 

32.020 

31.860 

31.760 

31.770 

31.850 

31.950 

32.030 

32.080 

32.120 

32.100 

32.060 

32.010 

31.960 



32.323 

32.392 

32.529 

32.618 

32.700 

32.789 

32.813 

32.803 

32.684 

32.485 

32.283 

32.155 

32.080 

32.042 

32.050 

32,109 

32.201 

32.287 

32.341 

32.369 

32.350 

32.301 

32.220 

32.145 



- 0.122 

-0.189 

-0.168 

- 0.120 

-0,129 

-0.133 

-0.153 

-0.134 

-0.075 

-0.093 

-0.135 

- 0.221 

-0.282 

-0.280 

-0.259 

-0.251 

-0.257 

-0.261 

-0.249 

-0.250 

-0.241 

- 0.210 

-0.185 



32,006 

32,065 

32.150 

32,239 

32,320 

32.412 

32.478 

32.470 

32.333 

32.182 

32.019 

31.845 

31.687 

31.594 

31,580 

31.621 

31.685 

31.780 

31.878 

31.941 

31,945 

31.933 

31.918 

31.902 



0.234 

0.205 

0,190 

0,211 

0.260 

0.248 

0.202 

0.180 

0.217 

0.228 

0.171 

0.175 

0.173 

0.166 

0.190 

0.229 

0.265 

0.250 

0.202 

0.179 

0.155 

0.127 

0.092 

0.058 



32.164 

32.228 

32.339 

32.428 

32.510 

32.600 

32.645 

32.636 

32.508 

32.333 

32,151 

32.000 

31.883 

31.818 

31.815 

31.865 

31.943 

32.033 

32.109 

32.155 

32.147 

31.117 

32,069 

32.023 



0.076 
0.042 
0.001 
0.021 
0.070 
0.060 
0.035 
0.013 
0.042 
0.076 
0.039 
0.020 
-0.024 
-0.058 
-0.045 
-0.015 
0.007 
-0.003 
-0.030 
-0.035 
-0.048 
-0.057 
-0 . 059 
-0.063 



RMS ERROR FOR FORWARD PREDICTION = 0.198 
RMS ERROR FOR BACKWARD PREDICTION = 0.198 
RMS ERROR FOR AVERAGE PREDICTION = 0,045 
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May 15, 1971 



ACTUAL 

DATA 

VALUE 

31.200 

31.200 

31.200 

31.190 

31.170 

31.170 

31.180 

31.230 

31.300 

31.360 

31.420 

31.430 

31.380 

31.280 

31.100 

30.920 

30.740 

30.590 

30.500 

30.490 

30.520 

30.640 

30.810 

30.950 



FORWARD 

PREDICTION 



VALUE 


ERROR 


31.350 


-0.150 


31.350 


-0.150 


31.360 


-0.160 


31.390 


-0.200 


31.450 


-0.280 


31.530 


-0.360 


31.600 


-0.420 


31.660 


-0.430 


31.740 


-0.440 


31,770 


-0.410 


31.710 


-0.290 


31.570 


-0.140 


31.410 


-0.030 


31.230 


0.050 


31.040 


0.060 


30.890 


0.030 


30.760 


-0,020 


30.670 


-0.080 


30.680 


-0.180 


30.740 


-0.250 


30.840 


-0.320 


30.980 


-0.340 


31.110 


-0.300 


31.210 


-0.260 



BACKWARD 

PREDICTION 



VALUE 


ERROR 


31.440 


-0.240 


31.570 


-0.370 


31.670 


-0.470 


31.720 


-0.530 


31.720 


-0.550 


31.670 


-0.500 


31,590 


-0.410 


31.550 


-0.320 


31.550 


-0.250 


31.570 


-0.210 


31,630 


-0.210 


31,690 


-0.260 


31.720 


-0.340 


31.730 


-0.450 


31.680 


-0.580 


31.550 


-0.630 


31.370 


-0.630 


31.230 


-0.640 


31.130 


-0.630 


31.020 


-0.530 


30.940 


-0.420 


30.940 


-0.300 


31.030 


-0.220 


31.180 


-0.230 



AVERAGE 

PREDICTION 



VALUE 


ERROR 


31.395 


-0.195 


31.460 


-0.260 


31.515 


-0.315 


31.555 


-0.365 


31.585 


-0.415 


31.600 


-0.430 


31.595 


-0.415 


31.605 


-0.375 


31.645 


-0.345 


31,670 


-0.310 


31.670 


-0.250 


31.630 


-0.200 


31.565 


-0.185 


31.480 


-0.200 


31.360 


-0.260 


31.220 


-0.300 


31.065 


-0.325 


30.950 


-0.360 


30.905 


-0.405 


30.880 


-0.390 


30.890 


-0.370 


30.960 


-0.320 


31.070 


-0.260 


31.195 


-0.245 



0.261 

0.440 

0.321 



RMS ERROR FOR FORWARD PREDICTION = 
RMS ERROR FOR BACKWARD PREDICTION = 
RMS ERROR FOR AVERAGE PREDICTION = 
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May 24, 1971 



FORWARD 
PREDICTION 
VALUE ERROR 

-0.051 



BACKWARD 
PREDICTION 
VALUE ERROR 



AVERAGE 
PREDICTION 
VALUE ERROR 



ACTUAL 

DATA 

VALUE 

30.940 

30.930 

30.980 

31.060 

31.180 

31.300 

31.370 

31.380 

31.340 
31.230 

30.990 
30.750 
30.560 
30.410 
30.330 
30.310 

30.340 
30.460 
30.640 
30.780 
30.880 
30.950 

30.990 
31.000 



30.991 

31.019 

31.112 

31.243 

31.363 

31.446 

31.451 

31.391 
31.274 
31.090 
30.897 
30.671 
30.483 

30.391 
30.355 
30.358 
30.420 
30.544 
30.682 
30.804 
30.891 
30.912 
30.904 
30.874 



-0.089 

-0.132 

-0.183 

-0.183 

-0.146 

-0.081 

- 0.011 

0.066 

0.140 

0.093 

0.079 

0.077 

0.019 

-0.025 

-0.048 

-0,080 

-0.084 

-0.042 

-0.024 

- 0.011 

0.038 

0.086 

0.126 



31.013 

31.012 

31.028 

31.067 

31.154 

31.265 

31.367 

31.452 

31.430 

31.315 
31,140 
30,902 
30.669 
30.465 

30.316 
30.251 
30.235 
30.294 
30.416 
30.570 
30.687 
30.780 
30.833 
30.840 



-0.073 

-0.082 

-0.048 

-0.007 

0.026 

0,035 

0.003 

-0.072 

-0.090 

-0.085 

-0.150 

-0.152 

-0.109 

-0.055 

0.014 

0.059 

0,105 

0,166 

0.224 

0.210 

0.193 

0.170 

0.157 

0.160 



31,002 

31.015 

31.070 

31,155 

31.258 

31.355 

31.409 

31.421 

31.352 

31.202 

31.018 

30,786 

30.576 

30.428 

30,335 

30.304 

30.327 

30.419 

30.549 

30.687 

30.789 

30.846 

30.868 

30.857 



-0.062 

-0.085 

-0,090 

-0.095 

-0.078 

-0.055 

-0.039 

-0.041 

- 0.012 

0.028 

-0.028 

-0.036 

-0.016 

-0.018 

-0.005 

0.005 

0.012 

0,041 

0.091 

0.093 

0.091 

0.104 

0.121 

0,143 



RMS ERROR FOR FORWARD PREDICTION = 0.094 
RMS ERROR FOR BACKWARD PREDICTION = 0.121 
RMS ERROR FOR AVERAGE PREDICTION = 0,070 
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May 30, 1971 



ACTUAL 

DATA 

VALUE 

30.300 

30.380 

30.440 

30.470 

30.470 
30.450 

30.420 

30.420 
30.440 
30.460 
30.530 
30.600 
30.630 
30.620 
30.540 
30.380 
30.210 
30.050 
29.930 
29.900 
29.920 
29.990 
30.080 
30.180 



FORWARD 

PREDICTION 



VALUE 


ERROR 


30.390 


-0.090 


30.380 


0.000 


30,340 


0.100 


30.320 


0.150 


30.320 


0.150 


30.340 


0.110 


30.390 


0.030 


30.470 


-0.050 


30.550 


-0.110 


30.610 


-0.150 


30.640 


-0.110 


30.620 


-0.020 


30.530 


0.100 


30.380 


0.240 


30.200 


0.340 


30.040 


0.340 


29.910 


0.300 


29.860 


0.190 


29.840 


0.090 


29.900 


0. 0 00 


30.010 


-0.090 


30.110 


-0.120 


30.190 


-0.110 


30.260 


-0.080 



BACKWARD 



PREDICTION 



VALUE 


ERROR 


30.380 


-0.080 


30.450 


-0.070 


30.490 


-0.050 


30,510 


-0.040 


30.510 


-0.040 


30.460 


-0.010 


30.350 


0,070 


30.250 


0.170 


30.200 


0.240 


30.190 


0.270 


30.220 


0.310 


30.260 


0.340 


30.330 


0.300 


30.380 


0.240 


30.420 


0.120 


30.450 


-0.070 


30.420 


-0.210 


30.340 


-0.290 


30.260 


-0.330 


30.170 


-0.270 


30.090 


-0.170 


30.080 


-0.090 


30.090 


-0.010 


30.140 


0.040 



AVERAGE 
PREDICTION 
VALUE ERROR 

-0.085 
-0.035 
0.025 
0.055 
0.055 
0.050 
0.050 
0.060 
0.065 
0.060 
0.100 
0.160 
0.200 
0.240 
0.230 
0.135 
0.045 
-0.050 
- 0.120 
-0.135 
-0.130 
-0.105 
-0.060 
- 0.020 



30.385 

30.415 

30.415 

30.415 

30.415 

30.400 

30.370 

30.360 

30.375 

30.400 

30.430 

30.440 

30.430 

30.380 

30.310 

30.245 

30.165 

30.100 

30.050 

30.035 

30.050 

30.095 

30.140 

30,200 



RMS ERROR FOR FORWARD PREDICTION = 0.158 
RMS ERROR FOR BACKWARD PREDICTION = 0.195 
RMS ERROR FOR AVERAGE PREDICTION = 0.113 
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